###### To do ######
# fails with data only argument + map (i.e., h or single color in link)
## parse map so it is possible to specify data directly in link()?
# add an easy way to generate sensible PDFs
# add option to include kernel density estimate with points (hitogram too?)
# More efficient stack algorithm
# improve output to integrate with other graphics packages

StackedDots <- function(data, map = link(), FUN = median,
  symm = FALSE, plot = TRUE, ...) {

  if (missing(map)) {
    map <- expression(link(x = data))
    data <- data.frame(data)
  }

  dat <- eval(expr = substitute(map), envir = data,
    enclos = parent.frame())

  stacks <- StackAlgorithm(data = dat@original[dat@ok],
    h = dat@h, interpolate = FUN)

  if (symm) {
    stacks$y <- with(stacks, ave(y, x, FUN = function(z) z - mean(z)))
  }

  dat@x <- stacks$x
  dat@y <- stacks$y

  stopifnot(ValidDotStackObject(dat))

  if (plot) {
    plot(dat, ...)
    return(invisible(dat))
  } else {
    return(dat)
  }
}

link <- function(x, h, color, shape = NULL) {
  key <- seq_along(x)
  ok <- key[!is.na(x)]
  xname <- deparse(substitute(x))

  if (missing(color)) {
    color <- "black"
  }

  if (missing(h)) {
    h <- length(ok)^(-1/2)
  }

  output <- new(Class = "DotStack",
    original = x,
    x = vector(mode = "numeric", length(ok)),
    y = vector(mode = "numeric", length(ok)),
    color = color, shape = shape,
    key = key, ok = ok, h = h, name = xname)

  return(output)
}

StackAlgorithm <- function(data, h, interpolate = median) {
  x <- matrix(FALSE, length(data), length(data),
    dimnames = list(seq_along(data), seq_along(data)))
  for (i in seq_along(data)) {
    x[, i] <- abs(data - data[i]) < (h * 2)
  }

  dat <- data.frame(stack = vector("character", length(data)),
    row.names = dimnames(x)[[2L]], stringsAsFactors = FALSE)

  ## this is computationally innefficient so I used
  ## .Internal calls directly here to speed it up
  while (length(x) > 0) {
    key <- .Internal(which.max(
      .Internal(colSums(x, dim(x)[1L], dim(x)[2L], TRUE))))
    neighbor <- .Internal(which(x[, key]))
    key.name <- dimnames(x)[[c(2L, key)]]
    neighbor.name <- dimnames(x)[[1L]][neighbor]
    dat[neighbor.name, "stack"] <- key.name
    x <- x[-neighbor, -neighbor, drop = FALSE]
  }

  dat$stack <- as.integer(dat$stack)

  cat(length(unique(dat$stack)),
    "iterations to make stacks", fill = TRUE)

  dat$x <- ave(x = data, dat$stack, FUN = interpolate)
  dat <- dat[order(dat$stack, data), ]
  dat$y <- sequence(rle(dat$stack)$lengths) * h * 2 - h
  return(dat[order(as.integer(rownames(dat))), c("x", "y")])
}
